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ABSTRACT 



Context. Weak gravitational lensing is one of the most promising tools to investigate the equation-of-state of dark energy. In order to 
obtain reliable parameter estimations for current and future experiments, a good theoretical understanding of dark matter clustering 
is essential. Of particular interest is the statistical precision to which weak lensing observables, such as cosmic shear correlation 
C^l functions, can be determined. 

Aims. We construct a fitting formula for the non-Gaussian part of the covariance of the lensing power spectrum. The Gaussian 
contribution to the covariance, which is proportional to the lensing power spectrum squared, and optionally shape noise can be 
^~ 5 included easily by adding their contributions. 

Methods. Starting from a canonical estimator for the dimensionless lensing power spectrum, we model first the covariance in the halo 
model approach including all four halo terms for one fiducial cosmology and then fit two polynomials to the expression found. On 
large scales, we use a first-order polynomial in the wave-numbers and dimensionless power spectra that goes asymptotically towards 
1.1 C p t for £ — > 0, i.e., the result for the non-Gaussian part of the covariance using tree-level perturbation theory. On the other hand, 
for small scales we employ a second-order polynomial in the dimensionless power spectra for the fit. 

Results. We obtain a fitting formula for the non-Gaussian contribution of the convergence power spectrum covariance that is accurate 
to 10% for the off-diagonal elements, and to 5% for the diagonal elements, in the range 50 < C < 5000 and can be used for single 
source redshifts z s e [0.5, 2.0] in WMAP5-like cosmologies. 

Key words, gravitational lensing - Methods: yV-body simulations - Cosmology: theory - large-scale structure of the Universe 
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Weak gravitational lensing by the large-scale structure, or cosmic shear, is an important tool to probe the mass distribution in the 
Universe and to estimate cosmological parameters. The constraints it provides are independent and complementary to those found 
by other cosmological probes such as cosmic microwave background (CMB) anisotropies, supernovae (SN) type la, baryon acoustic 
oscillations (BAO) or galaxy redshift surveys. The cosmic shear field quantifies the distortion of faint galaxy images that is induced 
by continuous light deflections caused by the large-scale structure in the Universe (e.g., Bartelmann & Schneider|2001||Schneider 
2006). Since this effect is too small to be measured for a single galaxy, large surveys with millions of galaxies are required to 
detect it in a statistical way. The cosmic shear signal has been successfully measured in various surveys, since the first detections 
of|Bacon efaL|(|2000); |Kaiser et ah] (|2000|>; |Van Waerbeke et al.| (|2000)>; [Wittman efaX] (|2000|i. Most recently, shear two-point 



correlation functions were measured in the Canada-France-Hawaii Te lescope Legacy Survey (CFHTLS) and were used to constrain 
the amplitude of dark matter clustering, cr 8 Q^ 6 , with 5% uncertainty (Fu et al. 2008} . 

The next generation of galaxy surveys will greatly improve the precision with which weak lensing effects can be measured 
( jAlbrecht et al.|2006] l enabling us to obtain, with accurate redshift information and tomographic measurements, precise constraints 
on the evolution of dark energy. However, the expected improvement in future data only leads to a significant improvement of 
the precision and accuracy of the cosmological interpretation if the systematic errors, the underlying physics, and the statistical 
precision of cosmic shear estimators are well understood. Systematics currently identified arise mainly from non-cosmological 
sources of shear correlations, i.e., intrinsic alignments of galaxies (e.g., |Schaefer 2008] for a review), and biases on the shear 
measurement (Mas sey et al.|2007 Sembolo ni et al.|20 08 ). This paper addresses the issue of the statistical precision of cosmic shear 
estimators, determined by the covariance of the estimator. Since much of the scales probed by cosmic shear lie in the non-linear 
regime, being affected by non-linear clustering, the covariance depends on non-Gaussian effects and has a non-Gaussian, as well as 
a Gaussian, contribution. Indeed, even though the non-Gaussianity of the shear field is weaker than that of the matter field due to the 
projection along the line-of-sight, various studies indicate that the non-Gaussian contribution to the covariance cannot be neglected 
when constraining cosmological parameters with weak lensing ( Scoccimarro et al.|199 9; Whit e & Hu|20 00; Coor ay & Hu|2 001 



Kilbinger & Schneid er]20Tj5] |Sembolon i et al.|2007||Takada & Jain|2009| l 

Most cosmic shear results are based on the measurement of two-point correlation functions of the shear field. Since, in general, 
the number of independent measurements is insufficient to infer the complete covariance directly from observations, one may 
derive it from ray-tracing maps of numerical A^-body simulations. This, however, requires a large number of realizations and, in 
addition, is very time-consuming if an exploration of the covariance in the parameter space is needed. An alternative is to compute 
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the covariance with an analytic approach. For shear two-point correlation functions, Schne ider et al.| ( |2002| ) derived an expression 
for the Gaussian contribution to the covariance. Semboloni et al. (2007 ) fitted the ratio between that expression and a covariance 
computed with A^-body simulations, containing both Gaussian and non-Gaussian contributions, providing thus a formula to compute 
the total covariance from the Gaussian term. In Fourier space, large-scale modes are independent and, differently from real space, 
the Gaussian contribution to the covariance of the convergence power spectrum (i.e., of the Fourier transform of the two-point shear 
correlation function) is diagonal and can be computed from the convergence power spectrum alone (Kaiser 1992, Joachimi et al 



2008), whereas the non-Gaussian contribution can be computed from the trispectrum of the convergence ( jScoccimarro et al.|1 999). 
The trispectrum, on large scales, can be accurately derived in tree-level perturbation theotyQ and, on small scales, is well represented 
by the one-halo term of a halo model approach. A non-Gaussian part of the covariance consisting of a perturbation theory term and 
a one-halo term was used, e.g., in Takada & Jain (2009). 

This paper aims at producing an accurate expression for the non-Gaussian contribution of the covariance of the convergence 
power spectrum that is fast to compute, contributing thus to accurate estimates of cosmological parameters. Following |Scoccimarro| 
|et al.| ( [T999] ) and Cooray ~& Hu| ( p0 01 ), we start from a canonical estimator of the dimensionless convergence power spectrum and 
use it to derive an analytic expression for the corresponding covariance. The various spectra involved are evaluated using the halo 
model approach of dark matter clustering ( |Seljak|2000"}|Ma & Fry|2000| [Scoccimarro et al.|2001||Cooray & Sheth|2002| ). The halo 
model approach assumes that all dark matter in the Universe is bound in spherical halos, and uses results from numerical A^-body 
simulations to characterize halo properties such as their profile, abundance and clustering behavior. 

The evaluation of the covariance of the convergence power spectrum in the halo model approach is time-consuming. In addition, 
it may be needed to repeat it for different cosmological models for the purpose of parameter estimation. To allow for a faster 
computation, we construct a fitting formula for the non-Gaussian part of the convergence power spectrum covariance. On small 
scales, we fit the halo model result with a polynomial in the non-linear dimensionless convergence power spectrum. On large scales, 
we fit the ratio between the halo model covariance and the perturbation theory covariance. We stress that it is a fit to the halo 
model covariance, not involving a covariance computed from A^-body simulations. The result is, however, calibrated by A^-body 
simulations, since they determine the halo model parameters. 

The paper is organized as follows. We define in Sect.[2]the reference cosmology, considering the growth of matter perturbations. 
We introduce the convergence spectra, construct an estimator for the dimensionless convergence power spectrum, and derive an 
expression for its covariance in Sect. [3] In Sect. [4] we describe the halo model approach, and compute the covariance of the power 
spectrum. The covariance depends on the values of halo model parameters, which are also defined here. It also depends on the 
power spectrum, bispectrum and trispectrum of the correlations of halo centers. Expressions for these spectra, in the framework 
of perturbation theory, are given in the Appendix. Section [5] tests the accuracy of the halo model predictions, for both the power 
spectrum and its covariance, against two sets of ray-tracing simulations. Section [6]presents the fitting formula for the non-Gaussian 
contribution to the covariance where its coefficients are given as function of source redshift. We conclude in Sect. [7] 



2. Structure formation in a ACDM cosmology 

Throughout this work we assume a spatially flat cold dark matter model with a cosmological constant (Q m + Q A = 1), as supported 
by the latest 5-year data release of WMAP results (Komatsu et al. 2009). The expansion rate of the Universe, H(a) = a /a, in such 
models is described by the Friedmann equation H 2 (a) = (fl m a~ 3 + Qa), where Hq = 100 h km s _I Mpc~' is the Hubble constant, 

Q m denotes the combined contributions from dark matter and baryons today in terms of the critical density p cr ; t = 3H^/(^nG), and 
Qa is the density parameter of the cosmological constant. The comoving distance to a source at a is then 

f 1 cda' 

where the scale factor is related to the redshift via the relation 1 + z = I /a using the convention a(to) = 1 today. 

In structure formation, the central quantity is the Fourier transform of the density contrast Six, t) = [p(x, t) - p(t)]/p(t), which 
describes the relative deviation of the local matter density p(x, t) to the comoving average density of the Universe p(t) at time t. We 
suppress the time dependence of 5 in the following. In this way, the mean density contrast is by definition zero, and we can describe 
matter perturbations in the early Universe as zero-mean Gaussian random fields. In this case, the statistical properties of the Fourier 
transformed density field, 

6(k) = j d 3 xe ikx 6(x), (2) 

are completely characterized by the power spectrum 

<§(*i)§(* 2 )> = (2nf5 D {k x + k 2 )P s {ki), (3) 

where (•) is the ensemble average and 6u denotes the Dirac delta distribution. Note that throughout this paper the tilde symbol is 
used to denote the Fourier transform of the corresponding quantity. 

In linear perturbation theory, which is valid on large scales, the power spectrum at a scale factor a is characterized by 

P lin (k, a) = Ak"> T 2 (k)D 2 (a) , (4) 



1 We refer by tree-level perturbation theory to the lowest, non-vanishing order of the considered quantity in perturbation theory. 
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where the amplitude A is normalized in terms of cr 8 , n s denotes the spectral index of the primordial power spectrum, and T(k) is the 
transfer function. Note that all Fourier modes of the matter density grow at the same rate, i.e., $(k, a) = 6(k)D(a), where 

H(a) r da' 

D(a) oc I (5) 

H Jo [a'H(a')/H ? 

is the growth factor which we normalize as D(a — 1) = 1. In the non-linear regime, i.e., on small scales, different Fourier modes 
couple and the Gaussian assumption cannot be maintained. Thus we have to consider higher-order moments of the density field to 
describe its statistical properties. In perturbation theory, it is possible to find analytic expressions for these moments, which hold up 
to the quasi-linear regime. In Appendix [A] we derive the expressions for the bispectrum and trispectrum in tree-level perturbation 
theory, which are the Fourier transforms of the three- and four-point-correlation functions, respectively. 



3. Covariance of the convergence power spectrum 



A central quantity in weak lensing applications is the two-dimensional projection of the density contrast 6(w6, w) on the sky, which 
is known as effective convergence k{6). It is obtained by projecting the density contrast along the backward-directed light-cone of 
the observer according to 



/"HI 
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Jo 



dw w G(w) 6(w6, w) , 



(6) 



where w = w(z) denotes the redshift-dependent comoving distance, wh is the comoving distance to the horizon and the weight 
function G(w) takes into account the distribution of source galaxies along the line-of-sight. We assume for simplicity that all 
background sources are situated at a single comoving distance w s = Vf(z s ), such that the weight function has the form 



G(w) 



■fin 



W s - W 



H(w s - w) . 



(7) 



where H(x) denotes the Heaviside step function. To take advantage of the Fourier properties, we analyze the statistical properties of 
the Fourier counterpart of k(6). For the theoretical consideration of the convergence power spectrum covariance we need the second- 
and the fourth-order moments, as will become apparent later. They are defined by 



(8) 
(9) 



where the subscript 'c' refers to the connected part of the corresponding moment and . ; = + . . . + 1 j is a sum of Fourier wave- 
vectors. The convergence power spectrum and trispectrum are calculated using the flat-sky and Limber's approximation ( |Kaiser| 
[T998)|Scoccimarro et al.|1999)|Bernardeau et al.|2002) : 







dwG 2 (w)P s 



, w 



w 



Jo w z \ 



t\ {2 tl {4 



w w w w 



(10) 
(11) 



where Pg and Tg are the corresponding three-dimensional matter power spectrum and trispectrum (Fourier transform of the four- 
point correlation function). 

We are interested in estimating the dimensionless convergence power spectrum 



P l( (€) = (€ 2 /27f)P K (€) l 



(12) 



and the corresponding covariance for wave-vectors of different length i. A natural choice for the estimator of the dimensionless 
convergence power spectrum is ( |Scoccimarro et al.|1999| Cooray & Hu|2 001 1 Takada & Bridle 2007 1 



1 r d 2 £ I 1 

PJQ = K(t)K{-t) , 

Aj meCi AAA)2n 



(13) 



which is unbiased in the limit of infinitesimal small bin sizes, since (p K {£d) = Pk(&)- Here, A = 4n f s k y denotes the solid angle of a 
survey with a fractional sky coverage of / s k y , and the integration is performed over the Fourier modes lying in the annulus defined 
by (t — A(i/2 < € < li + A€(/2, where A(i is the width of the z'-th bin. We denote the integration area formed by the annulus as A r (^,)- 
The evaluation of the covariance of the estimator, Eq. ( |13) , results in an expression of the form 



(2nf_ 

AAA) 



2'P 2 K (e i )8t i t i + T K {{iJj) 



NG 



(14) 
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where < P K ((i) is given by Eq. (12i, and 5^. denotes the Kronecker delta. The second term in square brackets is the bin-averaged 
convergence trispectrum, 

i 2 a r d% t\i\ (1S) 



J .w Am Jm*. Am w H 



where T K is the non-linear convergence trispectrum as defined in Eq. ( 1 1 1. To derive Eq. ( 14 1, we made use of the definitions of the 
convergence power spectrum and trispectrum in Eqs. ^ and and of the discrete limit of the delta distribution (5d(0) — > A/(2n) 2 . 
The derived expression consists of two terms: a Gaussian part C G , which scales as th e convergence power spectrum squared and 
only contributes to the diagonal of the covariance matrix (see, e.g. Joachimi et al. 2008 i, and a non-Gaussian part C NG , which scales 



as the dimensionless bin-averaged convergence trispectrum and introduces correlations between the wave-vectors of different bins 
(see, e.g. Scoccima rro et aL"1|1999| [Cooray & Hu||2001[ ). Both terms are inversely proportional to the survey area A, but have a 
different behavior with respect to the bin width A{j. While the Gaussian term decreases with increasing bin size, the non-Gaussian 
term is independent of the binning, since the bin area cancels out after the integration. 

We can analytically perform one of the integrations of the bin-averaged trispectrum in Eq. (JT3J. First, note that it only depends 
on the parallelogram configuration of the convergence trispectrum, i.e., setting (% = —t\, (3 = (% and {4 = -I2 in Eq. ( fTT) . Also, 
if we choose an appropriate coordinate system for the integration over the wave-vectors, the problem becomes symmetric under 
rotations and we can parametrize the convergence trispectrum by the length of the two sides of the parallelogram l\ and £2 and the 
angle between them, cos <p = (t\ ■ tj)jl \t%. Hence, we define 

TMJi,co&<p) = T K {t u -t u t 2 ,-h) . (16) 

Making use of the symmetry properties of this problem, one angular integration becomes trivial and the integration in Eq. (15i 
simplifies to 

T K (( i ,€ j )=— f —^-(\ f t-^A f d^T^i.^.cos^). (17) 
2?r J[f,|eft A T (ii) Jihieij A T ({j) Jo 

If the bin-width tS.ii is sufficiently small (A^; <c the integration area is A r (£) = 2n£tA€i, and we can make use of the mean value 
theorem. In this way, we approximate the integral in Eq. ( 17 1 by an angular average 

t 2 P. 



WiJj)^^- f d v -f^r,(4^,cos^). 

In Jo (2tt) 2 



(18) 



Note that if T K (€i,€j, cos (p) is independent of the angle between t; and £j, an approximation of the covariance can be calculated 
without having to perform an integration at all. In particular, this is the case for the 1-halo term of the three-dimensional matter 



trispectrum, as we will see later in Eq. (42 1 



4. Halo Model 



We have seen that the covariance of the dimensionless convergence power spectrum estimator consists of two terms: a Gaussian part, 
which is proportional to the dimensionless convergence power spectrum squared and a non-Gaussian part, which is the bin-averaged 
dimensionless convergence trispectrum (see Eq. 14 1. We will compute these terms using the halo model approach (Seljak 2000 Ma 
|& Fr y 2000; Scocc imarro et al.|2001| see also the comprehensive review by |Cooray & Sheth 2002). 



4.1. Overview 

With the assumption that all dark matter is bound in spherically-symmetric, virialized halos, the halo model provides a way to 
calculate the three-dimensional polyspectra of dark matter in the non-linear regime. In Sect. |4.5| below, we summarize the equations 
one obtains for the dark matter power spectrum and trispectrum. 

In the halo model description, the density field at an arbitrary position x in space is given as a superposition of all halo density 
profiles such that 

N N 

p(x) = 2^f(x- xr, nii, a) = 2^ u(x - xf, nii, a) , (19) 



where f(x—Xi\ nii, c,) denotes the density profile of the z'-th halo with center of mass at and u = f '/nii is the normalized profile. By 
parametrizing the halo profile in this way, we assume that the shape of the z'-th halo depends only on the halo mass m, and the halo 
concentration parameter c,-, which we define below. The dark matter polyspectra of the density field p(x) follow then from taking 
the ensemble averages, (X), of products of the density at different points in space. Assuming that the number of halos is N = hV, 
where n is the average number density of halos and V the considered volume, we compute the ensemble averages by integrating 
over the joint probability density function (PDF) for the halos that form the field ( Smith & Watts 2005 ), i.e., 



(X) 



r\ N 



d Xj dm, dc, 



P(x\ 



.,Xff,mi, . . . ,niN, c\, . , . , cn)X , 



(20) 
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where p(-, •, •) denotes the PDF. If one considers that position and mass of a single halo are independent random variables, the PDF 
factorizes as 

1 n(m) 

p(x;m,c) — p{x)p{m)p{c\m) — — p(c\m) , (21) 

V n 

where n(m) is the halo mass function and p(c\m) is the concentration probability distribution for halos given a mass m. 
4.2. Ingredients 

The halo model approach provides a scale-dependent description of the statistical properties of the large-scale structure. On small 
scales, the correlation of dark matter is governed by the mass profiles of the halos, whereas on large scales the clustering between 
different halos determines the nature of the correlation. As there are a multitude of models to describe the behavior on different 
scales, and an even larger number of parameters one has to set judiciously, there exists no such thing as a unique halo model. In 
order to have reproducible results, it is therefore necessary to specify ones choice of parameters. For this work, we will adopt the 
following parameters for the halo model: 

1. The average mass of a halo is defined as the mass within a sphere of virial radius r v ; r as m = (47r/3)r 3 ir A v j r p, where A v ; r denotes 
the overdensity of the virialized halo with respect to the average comoving mass density p in the Universe. Typically, values 
for A v i, are derived in the framework of the non-linear spherical collapse model (e.g. Gunn & Gott|1972" l. Expressions valid for 



different cosmologies are summarized in Nakamura & Suto ( 1997 1. In our implementation, we use the results which are valid 
for a flat ACDM-Universe, i.e., 

A vil -(z) = 18tt(1 + 0.4093* 2 ' 71572 ) , (22) 
where x = (Q m ' - l) I/3 /(l + z). We find for our fiducial WMAP5-like cosmology A vir (z = 0) = 349. 

Af-body simulations suggest that the density profile of a halo follows a universal function. We choose to use the NFW profile 
fNavarro et al.|1997] ), which is in good agreement with numerical results and has an analytical Fourier transform. It is given by 

p(r, m) — — = , (23) 

(r/r s )(l + r/r s ) 2 ' 

where p s is the amplitude of the density profile and r s characterizes the scale at which the slope of the density profile changes. 
For small scales (r < r s ) the profile scales with par 1 , whereas for large scales it behaves as p oc r~ 3 . The Fourier transform of 
the NFW profile is 



u(k;m,c) = ^ d^xpix; m, c) e lkx j J* d 3 xp(x; m, c) = J" drAnr 



2 sin(fcr) p(r; m, c) 
kr m 



ln(l+c)- ' 



1 +c 



sin 77 [Si([l + c]tj) - Si(77)] + cos tj [Ci([l + c]tj) - Ci(rj)] - -^^4 , (24) 

(l+c)rjj 

where 77 = kr YU /c, we truncated the integration at r wn in the second step, and introduced the concentration parameter c = r v i T /r s 
in the third step. Additionally, we use for the sine- and cosine integrals the definitions 



f 00 cosf , s f* sint , 
Ci(x) = - df , Si(jc) = dt . 

Jx t Jo 1 



(25) 

Jo 1 

The abundance of halos of mass mat a redshift z is given by 



p din v 
m 2 dlnm 

where we introduced the dimensionless variable v = v(m, z), 



n(m, z) - -jTZZi y /( y ) ' ( 26 ) 



v(m,z) = jr— , (27) 

where D(z) denotes the redshift-dependent growth factor, <r 2 (m) is the smoothed variance of the density contrast, and (5 sc (z) 
denotes the value of a spherical overdensity that collapses at a redshift z as calculated from linear perturbation theory. In our 
work, we use the expression from Nakamura & Suto ( |1997| >, which is valid for a ACDM Universe, 



Ssc(z) = 3(1 ^ )2 3 [l - 0.0123 ln(l + x 3 )] , (28) 

where x = (Q m ' - l) 1/3 /(l + z). The quantity has only a weak dependence on redshift and we find 5 sc (z = 0) = 1.675 for our 
fiducial model. 

The advantage of introducing v is that part of the mass function can be expressed by the multiplicity function v/(v), which has 
a universal shape, i.e., is independent of cosmological parameters and redshift. In this work, we employ the Sheth and Tormen 
mass function ( |Sheth & Tormen|1999| l 



vf(v) = A(p) [l + (qv 2 r p ) J- v exp (-qv 2 /2) , 



(29) 



which is an improvement over the original Press-Schechter formulation (Pres s & Schechter|1974| l. We use the parameter values 
p = 0.3, q - 0.707, and amplitude A(0.3) = 0.322, which follows from mass conservation. 
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4. The concentration parameter c = r v i r /r s characterizes the form of the halo profile. From A^-body simulations one finds that the 
average, c, depends on the halo mass ( Bullock et al.|200T i like 



c(m, z) 



c* 

1 + z \m* 



(30) 



where m, = m*(z = 0) is the characteristic mass defined within the Press-Schechter formalism as S sc (z = 0) = <r(m*). In the 
following, we will use the values c„ = 10 and a = 0.2 as proposed by Takada & Jain (2003 ). This implies that more massive halos 



are less centrally concentrated than less massive ones. However, results from numerical A^-body simulations (Jing 2000; Bullock 
et al.|200T) indicate that there is a significant scatter in the concentration parameter for halos of the same mass. Furthermore, 
Jing| ( |2000| proposes that such a concentration distribution can be described by a log-normal distribution 



p(c\m)dc — 



1 



: exp 



(In c - In c) 



1(tI 



dine . 



(31) 



Typical values for the concentration dispersion range from cr\ KC - 0.18 to cr\ nc — 0.32 (Jing 2000) [Wechsler et al.|20 02). Note 
that the width of the distribution cr lnc is independent of the halo mass. The variation of the halo concentration can be attributed 
to the different merger histories of the halos (Wechsler et al. 2002} . We will analyze the impact of this effect on different 
spectra in Sect. |4.7| When we use only the mean concentration parameter, we have to replace the probability distribution of the 
concentration, needed for example in Eq. ( pi) , by a Dirac delta distribution 



p(c\m)dc = 6d(c - c) c d In c . 



(32) 



On large scales, the correlation of the dark matter density field is governed by the spatial distribution of halos. Since the clustering 
behavior of halos and matter density differ, one introduces the bias factors bfyn, z) such that 



6h(6) = 6h(x;m,z) = b\(m,z)S{x) ■ 



bj(m,z) o 
2 ' - S 2 (x) + ... 



(33) 



In this way, the halo density contrast, 6^(6), is expressed as a Taylor expansion of the matter density contrast, 6(x). The bias 
parameters are in general derived based on the Sheth-Tormen mass function introduced above. For the linear halo bias one 
obtains then 



b\(m,z) - 1 + 



q[v(m,z)Y - 1 



2p 



(34) 



^sc(z) 5Az)[l+(q[v(m,z)?yV 
where p and q match the values used in the mass function. Expressions for higher-order bias factors can be found, e.g., in 



Scoccimarro et al. (2001 1. Since they only have a small impact on the quantities employed here, we take into account only the 



first-order bias. In Fourier space we may then write 

S h (6) = 6h(k;m,z) = b\(m,z)5{k) . 



(35) 



To obtain the final correlation function, one has to perform integrations along the halo mass and optionally along the halo 
concentration, with limits formally extending from to oo. In practice, we use the mass limits m m ; n = 10 3 hr x M e and w max = 
10 16 h~ l M Q . Masses smaller than m m j n = 10 3 /T 1 M Q give no significant contribution to the considered quantities, while, due to 
the exponential cut-off in mass, masses larger than m max = 10 16 /i _1 M Q are rare. For the concentration, we employ the integration 
limits c m i„ = 1 and c max = 10 3 . 

Due to the cut-off in mass, the consistency relation ( Scoccimarro et al.|2001| l 



Ammn(m,z)b\(m,z) — 1 



(36) 



does not hold. To cure this problem we consider a rescaled linear bias such that b\{m,z) — > bi(m,z)/b norm (z), where b norm (z) is 
the result of the integral in Eq. ( |36] l. In this way, one ensures that the halo term with the largest contribution to the correlation 
equals the perturbation theory expression on large scales (see Fig. [TJ. 



4.3. Building Blocks 

Using the ingredients described in the previous section, it is possible to define building blocks, which simplify significantly the 
notation for expressing the polyspectra ( |Cooray & Hu|200"T| ): 

dm I dcn(m,z) p(c\m)\ — \ bj(m,z)[u(ki;m,c) ■ ■ ■ u(k:;m,c)] . (37) 



In the case i - 0, we additionally define bo = 1, for consistency. 



J. Pielorz et al.: A fitting formula for the lensing power spectrum covariance 



7 



4.4. Power spectrum 



We can now compute the power spectrum from Eq. (20 1. The result consists of two terms, the 1-halo and the 2-halo terms, P$(k) 
Pih(k) + P 2h (k) ( |Seljak|2000| >. They are given by 



l r"' m » , r Cm « 

P\h(k) — I 6mn(m)m I dc p(c\m) \u(k;m, c) 
P Jm min J Cmln 



Pihik) = 



i r 

P Jm n 



dmn(m)mbi(m) 



£ 



dc p(c\m) u(k; m, c) 



PUk) , 



(38) 



(39) 



where Pn n (k) denotes the linear perturbation theory power spectrum, defined in Eq. (|4]), and we use the ingredients summarized 
earlier-on. Note that, for convenience, we omit the redshift-dependence in the notation. Using the building blocks from Eq. ( |37"] >, 
these terms can be written in the following compact form : 



Pih(k) = M Q2 (k, k) , P 2h (k) = [Mi i (k)) 2 P lin (k) . 



(40) 



The 1-halo term, Pu,, denotes correlations in space between two points in the same halo, whereas the 2-halo term, P 2 h, takes into 
account correlations between two different halos. Hence, the 1-halo term is dominant on small scales and the 2-halo term is dominant 
on large scales. Note that the 2-halo term converges to the linear power spectrum on large scales because of the consistency relation 
of the first-order halo bias factor (see Eq. [36]> and of the limit u(k; m, c) — > 1 for k — > 0. 



4.5. Trispectrum 

We compute now the dark matter trispectrum in the halo model approach. As discussed in Sect. [3] only parallelogram configurations 
of the trispectrum wave-vectors contribute to the covariance of the convergence power spectrum. Restricting our calculations to these 
configurations, we obtain four different halo term contributions, 



T s (ki , k 2 , cos <p) = r ih + T 2h + T 3h + T 4h , 



(41) 



which are simpler than in the general case. In addition, we neglect terms involving higher-order halo bias factors since, on most 
scales, they provide only a small correction (Ma & Fry 2000, Takada & Jain 2003). We further note that the perturbative expansion 
of halo centers, used in the calculations, was shown to become inaccurate on non-linear scales ( |Smith et"aL 2007|l. The con tributions 
to the trispectrum take the following forms, using the compact notation of the building blocks (see Cooray & Sheth 2002 for the 
expression of the halo model trispectrum including higher-order bias factors): The 1-halo term, dominant on the smallest scales, is 



T lb = M 0A (k 1 ,-k 1 ,k 2 ,-k 2 ). 



(42) 



Til + r 22 , consisting of a term T^, which corresponds to correlations of three points 



The 2-halo term has two contributions, T 2b - , 

within one halo and a fourth point in a second halo, and a term 7|?, which describes correlations involving two points in a first halo 
and the other two points in the second halo. They read, 



T» = 2M n (k u k 2 ,k 2 )M n (k 1 )P lin (k 1 ) + 2M l$ {k u k u k 2 )M n (k 2 )PUk 2 ) , 
Til = M 2 n (k u k 2 )[P M (\ki + k 2 \) + P lm (\ki ~ k 2 \)] . 

The 3 -halo term is given by 

T ih = 2M l2 (ki,k 2 )Mn(ki)M n (k 2 ) [B pi (k u k 2 ,-h - k 2 ) + B vt (k u -k 2 , -*i + k 2 )] 
Finally, the 4-halo term, dominant on large scales, is 

T 4h = M 2 U )M 2 n (k 2 )T ?t (ki ,-k u k 2 ,-k 2 ), 



(43) 
(44) 

(45) 

(46) 



and describes correlations of points distributed in four different halos. Note that, like for the power spectrum, the 2-halo term is 
computed from the linear power spectrum. On the other hand, the 3- and 4-halo terms depend on B pt and T pt , respectively, which 
are the lowest-order, non-vanishing, perturbation theory contributions to the bispectrum and trispectrum. Both spectra are derived 



in Appendix A. 2 



4.6. Convergence spectra 



The convergence power spectrum and trispectrum, needed to evaluate the covariance in Eq. ( 14 1, are computed by project- 
>rding to Eqs. (10 1 and (Hi. Fig. [T] shows the dimensionless convergence trispectrum, defined as T"(£) = 
-() for a square configuration where all wave-numbers have a length {, where n € {lh, 2h, 3h,4h) denotes the 



{J, 



ing P s and T s , 
{ 2 /2n i/TJW 

contribution from the corresponding halo term. The plot shows the individual contributions to the dimensionless T" (the projections 
of Eqs. [42p6*| , illustrating on which scales the individual terms are important, as well as the total contribution of all halo terms 
T^ ot . We show, in addition, the dimensionless projected Tf, which closely follows the 4-halo term. We see that the commonly used 
approximation 7^°' ~ Tf + 7^ lh is accurate for large wave-numbers {I > 10 3 ) but has a deviation of about 20% from the complete 
trispectrum 7~ K tot for small wave-numbers ({ < 10 2 ). 
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Fig. 1. Square configuration of the dimensionless convergence trispectrum T"{€) against wave-number t for the range 10 < I < 10 
for n 6 {lh, 2h, 3h, 4h, tot, pt). The upper panel displays the four individual halo terms, the sum of all four trispectrum halo terms 
T^ ot (solid line) and the tree-level perturbation theory trispectrum Tf', as indicated in the key. The lower panel shows the ratio 
between the indicated contributions and the complete trispectrum. The double dashed line illustrates the corresponding ratio of the 
approximation T l ° l ~ Tf' + T^ h . Note that we consider the 4-halo term only in the upper panel since it resembles the term Tf on 
large scales. 



4.7. Stochastic halo concentration 



The previous results were computed using the deterministic concentration-mass relation of Eq. (30 1. We now analyze the impact of 



scatter in the halo concentration parameter c on the covariance of the convergence power spectrum, using the stochastic concentra- 



tion relation given by the the log-normal concentration distribution of Eq. (31 1. 

Cooray & Hu| ( |2001| l, analyzed the effect of a stochastic concentration on the three-dimensional power spectrum and trispec- 



trum and found that the behavior of the corresponding 1-halo terms were increasingly sensitive to the width of the concentration 
distribution for smaller wave-numbers k. Furthermore, the effect of a stochastic concentration relation was more pronounced for the 
trispectrum than for the power spectrum, since the tail of the concentration distribution is weighted more strongly in higher-order 
statistics. 

Performing the same analysis for the projected power spectrum and trispectrum, we find a similar trend as in the three- 
dimensional case, but with a smaller sensitivity to the concentration width of the distribution on small scales. For <j\ nc = 0.3, we 
find, in the case of the 1-halo term of the power spectrum, a deviation from a deterministic concentration relation of about 1% - 2% 
for wave-numbers larger than ( ~ 10 3 , whereas, for the 1-halo term of the trispectrum, the deviation is of the order of 10% - 15% 
in the same f-range. Thus, when considering the covariance of the convergence power spectrum, one should take into account the 
concentration dispersion in the 1-halo term of the trispectrum but can safely neglect it for the power spectrum. Additionally, we find 
that a stochastic concentration has only a small impact on the 2-halo terms of the power spectrum and trispectrum ( Pielor z]2008 1 



From this analysis, we expect the effect of a concentration distribution to be the strongest on the non-Gaussian part of the 
covariance, which depends on the trispectrum. To directly infer the impact of a concentration distribution on the covariance, we 
calculate the 1-halo contribution to the non-Gaussian covariance, i.e., we perform the bin averaging of Eq. ( |15) , 

for two different concentration dispersions cr lnc , using our halo model implementation. Fig.|2]shows the ratio 



'y;/: mt/ m , (48) 



C^({i,{j;o-\ nc ) 
C? h G (^;cr lnc = 0) 

where cr lnc = denotes the deterministic concentration relation, for two values cr inc = {0.15, 0.3} in two contour plots. In agreement 
with the previous results, the largest impact of the concentration dispersion occurs for wave-numbers larger than 1; ^ €j - 4000. For 
cr\ nc - 0.15, the deviation of the covariance C^ h G from the original deterministic concentration relation is small, of about 3% - 6%. 
The effect becomes non-negligible for <x lnc = 0.3, where we find a deviation of 12% (for lj ^ lj > 4000) to 25% (for (i ^ (j > 8000). 
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J ln c 
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Fig. 2. Contour plots of the ratio R(ij, ij) between the non-Gaussian covariance 1-halo term contribution (see Eq. 48 1 computed with 
a c 



a concentration dispersion cr lnc and with a deterministic concentration, against wave-numbers ({,, €j) ranging from = (j = 500 to 



8500. The left panel is for concentration dispersion cr lnc = 0.15, whereas in the right panel cr lne = 0.3. 



Table 1. Cosmological parameters used to set up the initial power spectrum, which determines how the simulation particles are 
distributed initially. The simulations employ either the BBKS ( Bardeen et al.|[l98 6) or the EH ( jEisenstein & Hu||1998| l transfer 
function. Convergence maps with source galaxies situated at z s = 1 and z s = 2 were produced for both sets of simulations. 



Simulation 


a m 


n A 


h 


n h 


0"8 


n s 


r 


Zs 


T(k) 


Virgo 


0.3 


0.7 


0.7 


0.0 


0.9 


1.0 


0.21 


1 (2) 


BBKS 


Gems 


0.25 


0.75 


0.7 


0.04 


0.78 


1.0 


0.14 


1 (2) 


EH 



Table 2. Parameters used for generating the two A^-body simulations considered in this paper and for producing the resulting 
convergence maps with the multiple-lens-plane ray-tracing algorithm (see e.g., |Jain et al.|2000| |Hilbert et al.||2008[ ). From left to 
right these are: the side length, L^ ox , of the cubic simulation box, the number of particles, N„ m , used for the simulation, their mass, 
m par , and the number of available convergence maps, N nYdp , with area A map . 



Simulation 


L hm /h-' Mpc 


AW 


w par / h-' M Q 


^map 


A map /(deg) 2 


Virgo 


141.3 


256 3 


1.4 X 10 10 


200 


0.25 


Gems 


150.0 


256 3 


1.4 x 10 10 


220 


16.00 



5. Comparison with yV-body simulations 



In S ect.|4] we computed the power spectrum and trispectmm of the dark matter fluctuations. Projecting them according to Eqs. ([TO » 
and (|1 and inserting the result in Eqs. \\2\ , ( [14) >, and ( 15 i, we obtained the covariance of the convergence power spectrum 
estimator. 

To test the accuracy of the halo model predictions for the statistics of the dark matter density field, we compare the dimensionless 
convergence power spectrum and the corresponding covariance, calculated in the halo model approach, with results from two 
different ray-tracing simulations. 



5.1. Virgo and Gems simulation 

For our comparison, we chose one simulation from Jenkins et al. ( 1998! and ten simulations from Hartlap et al. ( 2009| l, which 
we denote in the following as Virgo and Gems simulation, respectively. The Virgo simulation was carried out in 1997 by the 
Virgo-Consortium for a ACDM cosmology (see Tab. [TJl with N pm - = 256 3 particles in a periodic box of comoving side length 
L box = 141.3 /r 1 Mpc (see Tab. ||. It uses the PP-/PM-code HYDRA, which places subgrids of higher resolution in strongly 
clustered regions. Structures on scales larger than 2/ so f t * 40 hr x kpc can be considered as well resolved. The Gems simulations 
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were set up in cubic volumes of comoving side length Lb ox = 150 h~ l Mpc with 256 3 particles (see Tab. [2}. The cosmology chosen 
reflects the WMAP5 results (Komatsu et al. 2009) and t hus has a lower value for erg than the Virgo simulation. It uses the GADGET-2 
code to simulate the evolution of dark matter particles ( |Springel|2005| and has a softening length of 2/ so f t m 30 /T 1 kpc. 



5.2. Ray-tracing 



The output of numerical simulations are three-dimensional distributions of A^ pal particles in cubic boxes over a range of redshift 
values. In order to compare the results with the predicted convergence power spectrum from the halo model, we make use of the 
multiple-lens-plane ray-tracing algorithm (see e.g., Jain et al.|2000 Hilbert et al. 2008) to construct effective convergence maps. 
The basic idea is to introduce a series of lens planes perpendicular to the central line-of-sight of the observer's backward light 
cone. In this way, the matter distribution within the light cone is sliced and can be projected onto the corresponding lens plane. By 
computing the deflection of light rays and its derivatives at each lens plane, one simulates the photon trajectory from the observer 
to the source by keeping track of the distortions of ray bundles. In this way, the continuous deflection of light rays is approximated 
by a finite number of deflections at the lens planes. As a result, one obtains the Jacobian matrix for the lens mapping from source to 
observer and can construct convergence maps. 

For both simulations, a similar number of around 200 effective convergence maps were produced, with source galaxies situated 
at a single redshift of either z s — 1 or z s — 2 (see Tab. gj. The maps produced with the Gems simulation have an area of 16 deg 2 , 
while the ones from Virgo are much smaller, with 0.25 deg 2 . 



5.3. Convergence power spectrum 

To test the accuracy of the halo model approach in describing the non-linear evolution of dark matter, we compare the dimensionless 
projected power spectrum, computed in the halo model approach, to the ones estimated from the numerical iV-body simulations. 

The dimensionless convergence power spectra of the simulations are measured from the real-space two-dimensional convergence 
maps of length L map and grid-size A^in- For this, we first apply a Fast Fourier Transform^ to obtain k({) on each grid-point. Then, 
we estimate the power spectrum at a wave-number t for the &-th convergence map by averaging over all Fourier modes in the band 
{ b -{{\(- A{/2 <\e\<( + At/2], i.e., 

1 I 1 

p»\f) = V — \k(£)\ 2 , (49) 

where N v (() denotes the number of modes in each band ( b of width A{. Finally, the ensemble average of the dimensionless power 
spectrum is obtained by averaging over the results of the various convergence maps, i.e., P{€) = jf— 2^7 P^(ti). The error 

* V map K — 1 

bars for the power spectrum estimate are computed from the dispersion over the N mav convergence maps used, and are due to 
sample variance. Since modes with a length similar to the side length of the convergence map are only poorly represented, the 
sampling variance is largest for small t. The effect is stronger for the Virgo simulation than for the Gems simulation, since the Gems 
convergence maps have a much larger area A map . 

Figj3]shows a good agreement between the halo model and the A^-body simulation convergence power spectra. We also include, 
in Fig. BTthe convergence power spectra computed using the fitting formulae of [Peacock & Dod"ds ( 1996 1 and Smith et al. (2003 ). 



Both, the halo model and the fitting formulae, show a better agreement with simulations for the lower source redshift case. On 
small scales, the halo model and the Smith et al. fitting formula agree well with the simulations, whereas the Peacock-Dodds fitting 
formula has too little power, in particular in the case of the Gems simulation. On intermediate scales ( t 10 3 ), the halo model is less 
accurate than the Smith et al. fitting formula, suggesting that the halo model suffers from the halo exclusion problem on these scales, 



as described, e.g., in |Tinker et al. (2005 1. This means that, while in simulations halos are never separated by distances smaller than 
the sum of their virial radii, this is not accounted for in the framework of our halo model and is probably the cause for the deviation. 
The good agreement of the Smith et al. formula is not surprising, since it is based on simulations of similar resolution than the 
ones we consider here. In contrast, a similar comparison using the convergence power spectrum estimated with the Millennium Run 
( Hilbert et al.|2 008 ) clearly favors the halo model prediction over of the two fitting formulae, with both fitting formulae strongly 



underestimating the power on intermediate and small scales. 

The good overall accuracy of the halo model results was expected since its ingredients, such as the mass function and the halo 
profile, were obtained from A^-body simulations. We note that, in this analysis, we used a deterministic concentration parameter, 
since the use of a stochastic one has only a small effect on the small scales of the convergence power spectrum. 



5.4. Covariance of the convergence power spectrum 

The similarity between simulation and halo model power spectra was to some extent expected. The ability to make an accurate 
description of higher-order correlations provides a stronger test of the halo model. Due to its important role in calculating the error 
of the power spectrum and for parameter estimates, we focus in this section on the accuracy of the covariance of the dimensionless 
convergence power spectrum. 

We need an appropriate estimator for the power spectrum covariance of the simulations. As each simulation provides Af map 
different K-maps, we have N m - dp realizations of the power spectrum. From these we estimate the covariance by applying the unbiased 



2 For the FFT we use an algorithm from the GNU scientific library (see http : //www. gnu. org/software/gsl/ for more details). 



J. Pielorz et al.: A fitting formula for the lensing power spectrum covariance 



11 




Fig. 3. Dimensionless convergence power spectrum V K (€) against wave-number I for galaxy sources at redshift z s = 1 (upper panels) 
and z s = 2 (lower panels). Points with error bars show measurements from the two sets of numerical simulations: Gems (left plots) 
and Virgo (right plots). Both simulations use a similar particle setup but differ in the area of the /(-maps (see Tabs. [T] and |2j. They 
are compared with the corresponding results obtained with fitting formulae from Smith et al. (long-dashed line), Peacock-Dodds 
(short-dashed line), and the halo model predictions for a deterministic halo concentration (solid line). See text for a discussion. 

sample covariance estimator which has for our purpose the form: 



N - 1 

1 v map i 



k=\ Vma p 



U=i 



U=i 



(50) 



where Pi\it) is the projected power spectrum estimate of the k-th effective convergence map at a wave-number €j (see Eq.|49j). We 
evaluate C s ; m for both, Virgo and Gems simulations, for the case z s = 1 . 

The halo model results, Chaio, are calculated as described in Sect.[4]and include all terms of the non-Gaussian covariance. In their 
computation, we use our fiducial halo model with the ingredients summarized in Sect. 4.2 and use the cosmological parameters 
values corresponding to each simulation, given in Tab. [T] for the case z s = 1. We consider both, a deterministic concentration-mass 
relation and a stochastic one, with dispersion o-\ nc = 0.3 for the 1-halo term of the trispectrum. 

We compare the halo model with the simulations covariance matrices considering their relative deviation 



AC; A ChaloWi^j) C s i m (£i,{j) 



C, 



(51) 



We note that although the number of available convergence maps for our simulations (~ 200) is too small to obtain a percentage level 
accuracy for the covariance estimate (Takahashi et a l.|2009[ ), we assume that the covariance from the simulations is the reference 



12 



J. Pielorz et al.: A fitting formula for the lensing power spectrum covariance 



°ln c=° 



a |nc =0.3 



12000 



9000 



6000 



3000 




~i 1 1 r 



ft . 



1 

0.8 
0.6 
0.4 
0.2 


-0.2 
-0.4 
-0.6 
-0.8 

-1 



3000 6000 9000 12000 



3000 6000 9000 12000 



Fig. 4. Relative difference ACij/Cij (see Eq. Dill between the halo model prediction for the covariance of the convergence power 
spectrum and the results from the Virgo simulation (z s = 1) against wave-numbers (Sufj). The binning scheme is linear with a 
constant bin width of Al = 720 ranging from £q = 720 to £20 = 15120. The left panel displays the full covariance, including all four 
halo terms for the trispectrum, and uses a deterministic concentration-mass relation denoted by <xi nc = 0. The right panel illustrates 
the same covariance but with a stochastic concentration distribution of width <x lnc = 0.3 for the 1-halo term of the trispectrum. 




1500 3000 4500 6000 1500 3000 4500 6000 



li l| 

Fig. 5. Relative difference ACij/Cij (see Eq. pit between the halo model prediction for the covariance of the convergence power 
spectrum and the results from the Gems simulation (z s = 1) against wave-numbers (£i,£j). The binning scheme is linear with a 
constant bin width of Al = 90 ranging from £q = 90 to £20 = 7020. The left panel displays the full covariance, including all four 
halo terms for the trispectrum, and uses a deterministic concentration-mass relation denoted by <xi nc = 0. The right panel illustrates 
the same covariance but with a stochastic concentration distribution of width o"i nc = 0.3 for the 1-halo term of the trispectrum. 



one, and put it in the denominator of Eq. (51 1. The comparisons are displayed in Figs. [4] and [5] for the Virgo and Gems simulation, 



respectively. In the case of the Virgo simulation with cn nc = 0, the best agreement is a relative deviation of 20% - 40% found on 
intermediate scales between £ ^ 1500 and £ ^ 5000. For small scales {£ > 6000) the halo model underestimates the simulation 
by 60% or more. On large scales, i.e., for wave-numbers £ < 1500, the halo model overestimates the simulation by around 60%. 
However, since the simulation covariances suffer from a large sampling variance due to the small size of the convergence maps, the 
comparison is not meaningful on these scales. The agreement improves for <x lne = 0.3, in this case there is a larger range of scales 
where the deviation is small. For the Gems simulation the agreement is much better. There is now an interval of scales (£ < 1500) 
where the deviation is in the range 0% - 40%. Throughout all the scales probed, the deviation on the diagonal of the covariances is 
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Fig. 6. Ratio C pt /C NG against wave-numbers The covariance in tree-level perturbation theory resembles the non-Gaussian 

covariance on very large scales (i.e., (, = lj < 20) to approximately 80% - 90% except along the diagonal and the closest off- 
diagonals. Here a good approximation of the non-Gaussian covariance requires at least one additional halo term. For wave-numbers 
li = lj > 100, an accurate description of the non-Gaussian covariance requires at least the 1-halo term. 

Table 3. The fiducial model used for the fitting procedure. See text for details. 
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around 20% - 40%. However, for off-diagonal components, at small scales, the agreement is still poor. Finally, the improvement of 
considering cn nc = 0.3 is less significant than in the Virgo case. 

Although the halo model predictions strongly deviate from the simulation estimates of the covariance on small scales, this anal- 
ysis does not necessarily imply a poor accuracy of those predictions. Indeed, the simulation covariances estimated in this analysis 
have a strong scatter. In particular, in the case of the Gems simulation, we found from bootstrap subsamples of 50 convergence maps 
that the resulting covariances can deviate up to 20% from the average covariance of the complete sample (Pielorz 2008). This is in 
agreement with the results by Takahashi et al. (2009) who need to use 5000 simulations to obtain an estimate of the matter power 



spectrum covariance at a sub-percent level accuracy. 

There are however very recent indications that Eq. ( 14 1 indeed underestimates the covariance of the convergence power spectrum 
on small scales ( Sato et aTj |2009). In particular, sample variance in the number of halos in a finite field is not accounted for by 



Eq. ( 14 1. Indeed, the mass function yields a mean number density, but there are fluctuations on the number of halos due to the 
large-scale mass fluctuations. Sample variance in the number of clusters in a volume-limited survey (Hu & Kravtsov 2003) has 
been considered in cluster abundance studies (e.g. ,|Vikhlinin et al.|2009| ). In the halo model framework, the sample variance in the 
number of ha los was derived in Takada & Bridle] ( |2"007| ) and its contribution to the covariance of the convergence power spectrum 
was found, in Sato et al. (2009 1, to boost the non-Gaussian errors of a 25 deg 2 survey by one order of magnitude on scales I w 10 4 . 
The increase is reduced for larger survey areas. 



6. Fitting formula for the covariance of the convergence power spectrum 

Future weak lensing surveys will provide much more precise measurements of the convergence power spectrum. In order to obtain 
robust constraints on cosmological parameters, accurate estimates of both the power spectrum and its covariance are needed. 



6.1. Methodology 

The number of measured power spectra is, in general, insufficient to infer the complete covariance directly from observations. One 
has thus to derive it either from ray-tracing maps of numerical A^-body simulations or with an analytic approach. A drawback of 
the first method is that it requires a large number of realizations and, in addition, is very time-consuming if an exploration of the 
covariance in the parameter space is needed. 

In the previous sections, we derived the covariance, and in particular its non-Gaussian part, with an analytic approach. This 
computation is, however, time-consuming and it might be useful to obtain an accurate covariance in a faster way. A first approach 
would be to rely on stronger approximations. For example, a commonly used approximation consists on evaluating the non-Gaussian 
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Table 4. Best-fit values for the parameters (po, p\,pi) of the redshift-fit for the set of coefficients (aq, . . . , ag). The 3-parameter fit is 
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covariance from T K « T pt + T\^, instead of using the full trispectrum. We saw in Sect. |4]that this approximation recovers the full 
trispectrum for scales I > 10 3 but deviates by ~ 40% on scales t < 10 2 for square configurations (compare also with Fig.[6jl. 

An alternative approach that we consider in the following, is to find a fitting formula for the halo model covariance that can be 
subsequently used without the need for implementing the halo model. We will provide a fit only for the non-Gaussian part of the halo 
model covariance, since the Gaussian part only depends on the non-linear convergence power spectrum and can thus be accurately 
computed without relying on the halo model. The inclusion of non-Gaussian errors increases the total covariance and one might 
think of fitting the ratio between the non-Gaussian and Gaussian terms. However, this is not a good quantity to fit, since the Gaussian 
contribution is diagonal and binning-dependent. In contrast, a similar ratio was fitted in the real space, where the Gaussian term 



is non-diagonal and binning-independent, using a non-Gaussian contribution measured from yV-body simulations (Semboloni et al. 



2007 ). In Fourier space, there is some relevant analogous information contained in the tree-level perturbation theory trispectrum. 
Indeed, pursuing the analogy with the real-space fit, T pt is a non-diagonal and binning-independent quantity, with a lower amplitude 
than the full trispectrum, approaching it at large scales. We thus compute C pt , the covariance predicted in tree-level perturbation 



theory on large scales (see Eq. A. 27 and Appendix|A|for a detailed derivation) and compare it with our covariance computed with the 
full halo model trispectrum, C , defined in Eq. jl4) . The ratio between the two covariances is shown in Fig. [6] In agreement with 
Fig. [T] the covariance predicted by tree-level perturbation theory contributes only ~ 50% to the complete non-Gaussian covariance 
along the diagonal. On very large scales (I < 100), the ratio C NG /C pt lies between 1.1 and 2. On smaller scales, C pt decreases fast 
and is no longer useful for fitting purposes. 

This discussion motivates us to use two different fitting formulae to model the non-Gaussian covariance over the whole range of 
scales. On large scales (1 < i < 200), we model the ratio C NG /C pt as a polynomial in the wave-numbers, t\ , Cj, and the dimensionless 
non-linear convergence power spectrum More precisely, we assume 

C<(4 lj) = Cpttfi, (j) [1.1 + a { mm + "iU + a 2 'P K ({ min ) + a 3 T K ({ max )] , (52) 

where { m m = mm(ij, (j) and l m . dx = max(^,', £j) which ensures the symmetry property C<(£,, (j) — C<.{tj, £i). In this way, C< — > 
1.1 C pt for large scales as predicted by the halo model (see Fig. [51, whereas the non-linear clustering on smaller scales is encoded 
in the first-order polynomial in wave-numbers and power spectra. On small scales, i.e. for 200 < I < 8000, we model directly the 
non-Gaussian covariance amplitude using a second-order polynomial in the power spectrum, such that 

C>(£, (j) = a 4 ?Wmin) + fl 5 n(^max) + a b V\(t mm ) + flv^max) + ^(AnirWmax) • (53) 

To ensure a smooth transition from small to large wave-numbers, we consider a linear combination of the two matrices defined in 
Eqs. ( f52] > and ( |53~| with weightings of third-order in the wave-numbers, such that the full non-Gaussian covariance becomes 

C£ G (4 (j) = ; *° - - C<(€u ij) + J {i * {j \ ^ C^ti, £j) , (54) 



with the transition scale fixed at £q = 600. We use thus a model with nine free parameters, which is a compromise between 
expressive power and simplicity, and perform the fit using the range 1 < { < 8000. More precisely, we find the coefficients for the 



final fit formula by performing two least square fits to Eq. p2\ and Eq. (53 1, respectively. To obtain the coefficients for small (large) 
scales we computed the corresponding covariance as described above within the halo model approach with 61 bins in the range 
10 2 < I < 10 4 (1 < I < 10 3 ) and employed the condition € t + €j > 200 + €j < 200) for the fit. 

The covariances used in the fitting procedure, both tree-level perturbation theory and halo model covariance, depend on pertur- 
bation theory polyspectra. These were evaluated from the expressions derived in Appendix[A]using a linear matter power spectrum 
computed with the Eisenstein-Hu transfer function (Eisen stein & Hu|1998| l and assuming the WMAP5-like fiducial model shown in 
Tab. 3] The non-linear convergence power spectrum, used in the polynomial expressions, was evaluated from the same linear power 



spectrum. In addition, the halo model non-Gaussian covariance was evaluated using the input parameters as described in Sect. 4.2 
including a stochastic concentration-mass relation with cr lnc = 0.3 for the 1-halo term of the trispectrum. 
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Fig. 7. Coefficients a, of the fitting formula obtained for various source redshifts. Left (right) panel shows the coefficients of the 
large (small) scales fitting formula. For convenience, scaled versions of the coefficients are shown, as indicated in the key. Note for 
example that as, 07 and ag change sign. The lines show the redshift fit, Eq. ( f55) , with the parameters given in Tab. [4] They provide 
a good fit for z s > 0.5 (indicated by the vertical line). 



6.2. Redshift dependence 

All quantities, C NG , C pt , and V,,, were evaluated for several values of source redshifts (assuming a single source redshift plane), 
in the range 0.1 < z s < 2.0. For each redshift, we perform the fit and find a set of best-fit values for the coefficients (ao, . . . , as). 
Next, we fit the best-fit values of each coefficient as a function of redshift. Some of the best-fit values are increasing functions of 
the source redshift, while others decrease with redshift, and others still are non-monotonic. They all are, however, monotonic in the 
redshift range of interest for current and future weak lensing surveys, 0.5 < z s < 2.0. In this range, the best-fit values a,(z s ) are 
accurately fitted with 

a(z s ) = -^ +P2 . (55) 

Tab. [4] shows the resulting values of the 27 parameters, which define our fitting formula for the halo model covariance of the 
convergence power spectrum. The behavior of the coefficients with redshift is shown in Fig. [7] We note that a few of them change 
sign, with most remaining always positive. In addition, their amplitudes may be quite distinct, since they are applied to quantities of 
quite different amplitudes, such as wave-number, power spectrum or power spectrum squared. The coefficients of the fitting formula 
for large scales, Eq. ( |52) , all decrease with redshift in absolute values. This compensates the increase of the power spectrum with 
redshift, allowing for a decrease of the ratio C NG /C pt with redshift, as expected. 



6.3. Accuracy of the fitting formula 

To test the performance of the fitting, we calculate the relative deviation between the non-Gaussian covariance computed with the 
fitting formula and the halo model one. The deviation, 



AC iJ \_C™(e i Jj)-C™(£ i Jj) 



(56) 



is computed for every wave-number pair (fit, £j) and shown in Fig. [8] for the case of z s = 1- The upper left panel of Fig. [8] shows 
the absolute value of the deviation for each element of the covariance matrix, while the other panels show the deviation along cuts 
through diagonals of the covariance. 

The fit works quite well on the off-diagonal elements that are close to the diagonal, showing an average overestimation of 10%. 



When moving along any off-diagonal, from larger to smaller scales, (lower panels of Fig. (8), we move from the first fit, Eq. (52 1, 
where the deviation is mostly positive, to the second one, Eq. ((53), where the deviation is mostly negative. The transition occurs at 
the local maximum, which indicates that none of the fitting formulae should be extrapolated to the other region. The fits break down 
on the smallest scales, with the deviation increasing very rapidly when the largest scale reaches =s 5000. 

As we move away from the diagonal the fit gets increasingly worse, in particular the deviations are larger than 100% in the 
region shown in black in the upper left panel of Fig. [8] which correlates very small with very large scales. The reason for this is that 
this region was effectively not fitted, since it is not contained in any of the two blocks fitted by Eq. ( |54) . Restricting to the range 
where both scales are between 50 < I < 5000, roughly 90% of the elements show deviations between -25% and +25%, with the 
average of the absolute deviations being 10%. This range contains also 9% of outliers where the deviations are larger than +25%. 
The outliers occur in the low-amplitude correlations between the largest (50 < I < 200) and smallest (3000 < t < 5000) scales. 

The fit is worse on the diagonal than on the first off-diagonals (where lijlj < 50). On the diagonal, the fit always underestimates 
the covariance. Scales in the range 50 < i < 5000 show a deviation between -40% and -10%, with an average of -20%. The 
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Fig. 8. Accuracy of the fitting formula for z s = 1. Upper left panel: relative deviation ACy/Cy, defined in Eq. i56i. Other panels: 
diagonal cuts through the fitted non-Gaussian covariance, showing the covariance on the diagonal 1; = {j (upper-right panel), for 
li = 4.6^ (lower-left) and for €% = I0£j (lower-right). All three panels show the fitted and the halo model non-Gaussian covariances 
(as well as the Gaussian for the upper-right panel) as function of the lowest scale l> in the upper part, and the deviation in percent 
in the lower part, where the dashed lines mark the 25% level. For the Gaussian covariance we employed a logarithmic binning with 



6 1 bins in a range from £\ ovi — 1 to ( u 



10 3 



accuracy degrades at larger scales, which is not a problem since for I < 50 the non-Gaussian contribution to the total covariance is 
negligible, as seen in Fig.[8](upper right panel). Adding the Gaussian contribution to the fitted non-Gaussian one, the underestimation 
in the diagonal elements is always better than 10%, with an average of 5%. 

In summary, the fitting formula for the cosmic variance, including Gaussian and non-Gaussian contributions, has an average 
accuracy of 10% in the off-diagonal and 5% in the diagonal. It is valid when both scales are in the range 50 < t < 5000, corre- 
sponding to 2' < 6 < 5° in real space. This is roughly the range used in the latest results from CFHTLS-Wide ( |Fu et al.|20 08 ). This 
range includes the scales where non-Gaussianity is relevant, i.e., where the cosmic shear error budget is both dominated by cosmic 
variance and has important contributions from non-linear clustering (see Sect. |6.4|>. 



6.4. Impact of the accuracy of the fitting formula on parameter estimations 

We study the impact of the fitting formula accuracy on the estimation of cosmological parameters, using a Fisher matrix approach. 
For this, we need to take into account not only the Gaussian and non-Gaussian contributions to the covariance, but also the noise in 



the observed power spectrum. In practical applications, the convergence field in Eq. ( 13 i is obtained from the observed ellipticities of 
the source galaxies. The intrinsic ellipticity field (i.e., in the absence of a gravitational lensing effect) is assumed to have zero mean 
and rms of <x e per component. This shape noise contaminates the observed power spectrum. Assuming that the intrinsic ellipticities 
of different galaxies do not correlate, the shape noise contribution to the covariance of the power spectrum is diagonal and given by 
the new terms arising in Eq. ( fT~4| > when replacing the power spectrum, in that expression, by the observed one defined as ( |Kaiser| 
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Fig. 9. Left panel: Relative contributions to the diagonal of the convergence power spectrum covariance. The binning-dependent 
Gaussian to non-Gaussian ratio (solid line; see caption of Fig. [8]for the employed binning scheme) and the shape noise to cosmic 
variance ratio for two surveys with n = lOarcmin 2 and n = 40 arcmin -2 (dotted), are shown as function of multipole £. The ticks 
at the intersection of the 10% line with the ratio curves enclose the approximate range where the non-Gaussian contribution is 
important for the EucLiD-like survey. Right panel: Fisher ellipses in the (Q m , cr 8 ) plane (2<x contours) calculated with the halo model 
covariance (solid) and the fitted covariance (dashed) for the DES-like (large ellipses) and EucLiD-like (small ellipses) surveys. 
Compared to the halo covariance case, the fitted covariance error ellipses are decreased (enlarged) by 13% (10%) for the DES 
(Euclid) case. 



[1992) , 



(T 2 



pf\() = P K {£) + _1, (57) 



n 



where n is the number density of source galaxies. 

We consider the following three surveys: a medium-deep weak lensing survey covering an area of 170 deg 2 with n — 10 arcmin 2 , 
like the current CFHTLS-Wide, a wider survey of similar depth covering 5000 deg 2 with n = 10 arcmin , like the _planned Dark 
Energy Survey (DES^ and a wide and deep survey with 20 000 deg 2 and n = 40 arcmin 2 , like the proposed Euclid^ For all three 
we assume cr e = 0.3 and compute the covariance using the halo model and the developed fit formula in the range 50 < £ < 5000. 
Additionally, we add the covariance of a Gaussian contribution using a logarithmic binning with 61 bins in a range from flow = 1 to 
Aip = 10 s . The wider surveys will measure correlations on scales larger than £ = 50, which we do not consider here for comparison 
purposes. Note also that for very large scales the flat-sky approximation breaks down. 

The left panel of Fig. [9 compares the different terms contributing to the diagonal of the covariance, by showing the Gaussian 
to non-Gaussian ratio and the shape noise to cosmic variance ratio, where by cosmic variance we denote the sum of the Gaussian 
and non-Gaussian contributions. The ratio between the Gaussian and non-Gaussian terms is independent of the survey and, for 
our particular choice of binning, non-Gaussianity starts to affect the diagonal around £ = 300, where its amplitude is 10% of the 
Gaussian amplitude, and dominates from £ w 1000 onwards. Shape noise, including both pure shape noise and the coupling with 
cosmic variance, also becomes important on small scales, but in a survey-dependent way. It is as large as the cosmic variance on 
£ m 200 (£ ~ 1000) for surveys with 10 (40) galaxies per arcmin squared. The vertical lines in Fig. [9] (left panel) show, for the 
EucLiD-like survey, the range where the non-Gaussian contribution to the diagonal is non-negligible, i.e., where it accounts for more 
than 10% of the cosmic variance while having an amplitude of at least 10% of the shape noise. This range is roughly 300 < £ < 3000, 
or approximately 4' < 6 < 50'. This is a rough estimate of the minimum range where the fitting formula is required to have a good 
accuracy. 

In addition, the accuracy of the off-diagonal terms is crucial, since the non-Gaussianity is the sole contribution there. To evaluate 
the required range of validity of the fitting formula, in a way that includes the off-diagonal elements and is independent of bin width, 
we define the signal-to-noise ratio ( |Takada & Jain|2009 1, 



(|) =Z^(A)Cy 1 ^). (58) 

U 

For each survey, we compute the signal-to-noise ratio (SNR) using all scales between £ = 50, the largest scale where the fitting 
is valid, to successive values of ^ max . The SNR increases with £ max , as more scales are included in Eq. (58 I. The increasing rate, 
however, decreases with increasing f max , tending to zero when additional scales do not carry additional cosmological information. 
The value of f max for which this saturation occurs increases with decreasing shape noise, and gives a good indication of the range 
where the fitting formula is required to have a good accuracy in order to get accurate estimates of cosmological parameters. We find 

3 http : //www . darkenergysurvey . org/ 
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that, when increasing £ max from £ mdx = 5000 to £ max — 10000, the SNR increases by a factor of 1.03 for the CFHTLS-like survey 
and by a factor of 1.09 for the EucLiD-like survey. Hence, even for the latter there is no much gain in reaching t > 5000. 

We consider now the Fisher information matrix, which to first order, neglecting the cosmology dependence of the covariance 
matrix, is given by 

Cjl * 1 , (59) 



F a p - Yj 



dp a 11 dpp 



where the derivatives are taken w.r.t. a set of cosmological parameters p a . The Fisher matrix defines the error ellipsoid in parameter 

space, with \F~^ a ) yielding the lcr marginalized error on p a , and {l/F aa ) 1 ^ 2 giving the lcr error on p a assuming the other 
parameters are perfectly known. 

We compute Eq. ( |59| for both the fit and halo model covariances, for each of the three surveys. We perform the derivatives at 
the fiducial model of Tab. [3] varying only 2 cosmological parameters (Q m , cr%). For each survey, we compare the areas of the two 2<x 
error ellipses (which defines the inverse of the figure-of-merit) thus obtained. For all three cases, there is a good agreement between 
the two ellipses. 

For the two cases with large shape noise, CFHTLS and DES, we find that the fitting formula underestimates the Fisher ellipse, 
as compared to the halo model covariance. This is expected because with an enhanced diagonal the correlations between bins are 
weaker, and the result is dominated by the accuracy of the diagonal, where the fitting formula underestimates the covariance, as we 
saw earlier on. The deviation is however weak, the areas of the ellipses obtained using the fitting formula are 13% smaller than the 
halo model result, and the deviation is uniformly distributed on the parameter space (see Fig. |9j right panel, large ellipses). This 
implies that the deviation on the marginalized constraints is much smaller, and we find that the fitting formula underestimates the 
errors on both Q. m and cr 8 by only 1%, for both surveys. This corresponds to a deviation of 0.1% (0.01%) of the parameters values, 
for CFHTLS (DES). In contrast, for the Euclid case, where the covariance has larger correlations, the fitted covariance produced an 
ellipse slightly larger than the halo model one, by about 10% (see Fig. [9] right panel, small ellipses). 

6.5. Covariance of real -space estimators 

For practical purposes it is sometimes more convenient to study real-space correlations rather than correlations in Fourier space. We 
therefore define an estimator of a general second-order cosmic shear measure which is related to the convergence power spectrum 
estimator by 

f°° dtt 

m = w(W)PA{), (60) 

Jo 2tt 

where W(x) is an arbitrary weight function. A well-known example of this equation are the shear two-point correlation functions 
£+(0) and £_(0) with weight functions W{€0) = J (W) and W((6) = J 4 (W), respectively (see discussion in |Joachimi et al.||2008) . 
Using the definition Eq. ( |60| ), we find for the relation between the covariance of the real-space estimators to the covariance of the 
Fourier-space estimators: 

X°° dt £ r°° df £' 

~2n~ W{M) I ^W'^Covf^XP^')] , (61) 

which is related to the covariance of the dimensionless power spectrum used in the fitting formula by (2n) 2 Cov = 
Inserting the result of the dimensionless power spectrum covariance given in Eq. t fB] ) yields 

A„ r*°° Af> 1 Ap /"»°° AP' 

Cov [f (0), f (00] = — J o ^ P 2 K (()W(W)W(W') + - J o — W(ty J o — T K (£, C)W(i'6') . (62) 

We find that the Gaussian part of the real-space covariance is independent of the binning scheme and is non-diagonal in contrast to 
the covariance in Fourier space. 



7. Conclusions 



We present a fitting formula for the halo model prediction of the non-Gaussian contribution to the covariance of the dimensionless 
power spectrum of the weak lensing convergence. The formula was constructed assuming a ACDM cosmology with WMAP5-like 
cosmological parameters. In particular, it was obtained for Q m = 0.28, cr% - 0.82 and other parameter values as shown in Tab. [3] 
It is valid for a scale range of 50 < t < 5000, corresponding to 2' < < 5° in real space and can be used for surveys with galaxy 
source redshifts z s e [0.5,2]. In this range, it reproduces the results of a full implementation of the halo model approach, with a 
scale-averaged accuracy of 10% in the off-diagonal, and 5% in the diagonal elements. The formula also allows us to recover the 
halo model (£2 m , cr 8 ) error ellipses within 15%. The range of validity of the formula and its level of accuracy render it applicable to 
low shape noise scenarios from next generation weak lensing surveys. 
To use the formula, shown in Eq. ( 54 1, one needs three quantities : 



- The non-Gaussian contribution to the covariance of the convergence power spectrum in tree-level perturbation the ory C pt , shown 
in Eq. (A. 27 1. This involves the computation of the convergence trispectrum in tree-level perturbation theory, Eq. (A. 26 1, which 
requires the calculation of the linear power spectrum and of the Fj and Fj, coupling functions (Eqs . | A. 1 0| and [ATI 2\ . 
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- The non-linear convergence power spectrum. 

- The 9 coefficients of the fit, which are obtained by inserting the 27 values given in Tab. |4]in Eq. ( 55 I, for the required redshift. 



The Gaussian contribution, calculated from the non-linear convergence power spectrum as given in Eq. ( 14 1, may then be added to 
the result of the formula, to obtain the total covariance. This is the covariance of the estimator without noise, or the cosmic variance. 

The work presented in this paper is based on the assumption that the halo model is a powerful approach to probe non-linear 
clustering. We tried to test this assumption against results from N-body simulations, but our comparisons were inconclusive. Indeed, 
such analysis requires ray-tracing simulations with both large number of convergence maps and large convergence map area. Only 
then it would be possible to minimize the effect of sampling variance in the simulations, which hindered our attempted tests. Such 
simulations would also allow us to consider error bars for the estimate of the simulation covariance. There are however indications 
that the halo model approach underestimates the non-Gaussianity of the covariance and there are attempts to include additional 
contributions ( |Takada & Jain|2009[ |Sato et al.|2QQ9] > . 

The reliability of the halo model for higher-order polyspectra also needs to be studied in more detail. Some work in this direction 



are the analyses of the impact of the triaxiality of the halo profiles (Smith & Watts 2005; Smith et al. 2006), and halo exclusion 
effects ( |Tinker et al.|2005} . Moreover, the issue of halo substructure (Dolney ~et al.|2004] l and of the effect of a stochastic concen- 
tration parameter have to be understood properly. This paper also addresses this last issue. We analyze the impact of a stochastic 
concentration parameter on the covariance of the convergence power spectrum. We found that the effect can safely be neglected 
for the Gaussian contribution, with the convergence power spectrum varying only slightly, and at small scales, for concentration 
scatters of cr lnc =s 0.2 - 0.3. For the non-Gaussian contribution the effect is more pronounced due to the higher sensitivity of the 
trispectrum to a stochastic concentration relation. In the case of the 1-halo term of the trispectrum, we find it useful to take into 
account a concentration dispersion of cr\ ac > 0.3. The deviation to a deterministic concentration-mass relation is larger than 12% 
for wave-numbers ( > 3000. 

Although the fitting formula we obtained provides a more thorough estimate for the non-Gaussian contribution to the power 



spectrum covariance than the earlier approximation of Semboloni et al. d2007 ) (global accuracy of ~ 20% along the diagonal, which 



becomes less for the off-diagonals), there is still room for improvements. One drawback of our approach is that it requires the 
computation of the convergence trispectrum in tree-level perturbation theory. Furthermore, it is only applicable to a small range of 
WMAP5-like cosmologies and is only valid in the interval 50 < i < 5000. A possible way to avoid these problems and extend 
the accuracy of the fitting formula might be to construct it entirely from its three-dimensional counterpart, the three-dimensional 
covariance of the matter power spectrum. This has the advantage that perturbations of different length scales are not additionally 
mixed due to projection effects, which might allow us to cover a wider range of cosmologies. The desired projected covariance 
could then be obtained by performing an additional integration along the redshift-space. We will address this issue in a future paper. 
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Appendix A: Cosmological Perturbation Theory 

On large scales, different Fourier modes evolve independently from each other and thus conserve the Gaussian behavior of the 
density perturbation field 6(k, a). It is therefore convenient to work in Fourier space and Fourier transform the fields as well as 
the non-linear fluid equations (consisting of continuity, Euler and Poisson equation) that describe their evolution in an expanding 
Universe. Contrary to linear perturbation theory, there is a coupling between different Fourier modes mediated by the coupling 
function a(k\,k 2 ) on smaller scales. In this case, the fluid equations for the density contrast and the irrotational peculiar velocity 
field 6 — V ■ u in Fourier space are given by (e.g., Bernarde au et al.|2002 1 



J (2tt) 3 J 



aS(k,a) + e(k,a) = - <rk 2 6 D (k - k x - k 2 )a(k u k 2 )0(k u a)6(k 2 ,a) , (A.l) 



2a J (2n) 

where we introduced the two fundamental mode coupling functions 



a6(k,a)+ ""»"'" l(k,a) = - f f d 3 k 2 S D (k - ki - k 2 )B{k u k 2 )8{kud)6{k 2 ,a) , (A.2) 

2a J (2/r.r J 



(k l + k 2 )-k l \k\ + k 2 \ z (ki ■ k 2 ) 
a(k u k 2 )= - 2 , P{ki,k 2 )= — -— 2 . (A.3) 

For an Einstein-de Sitter (EdS) cosmology it is possible to find a perturbative ansatz that separates the scale- and time dependencies, 
whereas for a general ACDM model it is impossible to find a separable solution to Eqs. ( |A.1[ ) and ( |A.2[ ). However, |Scoccimarro et aL\ 
( 1998 ) showed that it is possible to find a separable solution in any order if one makes an approximation that is valid at percentage 
level. One indeed finds then the same recursion relation as in the EdS case. The ansatz is then 

oo oo 

6(k, a) = Y D n (a)6 n (k) , 6(k, a) = -aV D n (a)0 n (k) . (A.4) 

/?= 1 n=\ 

Therefore, we find that the whole information on cosmological parameters is encoded in the growth function due to its dependence 
on the Hubble parameter (see Eq.|5]). 
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A.1. Coupling Functions 

The n-th order density contrast and the divergence of the peculiar velocity in Eq. (|A.4|i is given by 



/d q\ d q„-\ f o - 
n^'"~n^f J ^^dC*-?!...™)^^!,...,^!^)---^^), (A.5) 

§*(*) = J |^ • • ■ J d 3 q n 6 D (k - q x ... n )G n {q„- ■ • , ft)§i(ft) ■ • ■ §i(ft) , (A.6) 

and the n-th order coupling functions F„ and G„ are obtained by the following recursion relations ( Jain & Bertschinger|1994| l 



H-l 



F n (q x ,. ..,q n ) = Y j [(2n + l)ar(*i, k 2 )F n . m (q m+l , . ..,q n ) + 2fi(k u k 2 )G n . m (q m+u . . . , q„)] , (A.7) 

i (2n + 3)(n - 1) 



m=l 
n-1 



G n (q u . .. , ft) = Y S ! ^P i fW*i i k 2 )F n - m (q m+l ,. ..,q„) + 2np{k u k 2 )G n . m (q m+1 , . . . , q„)] , (A.8) 
(2« + 3)(n - 1) 

where £i = q 1 + . . . + q m and £ 2 = ft,+i + ■ • • + q„- The initial conditions for these recursion relations are F\ = 1 and G\ = 1. To 
get the functions F ( „ k) and G„ s) that are symmetric in its arguments, one must perform the following symmetrizing procedure 

f<%!,. ..,?„) = — F n(q„ m , ■ ■ • , q„(n)) > c^Cft, • ■ ■ . q n ) = ^ X G " ( ^a>' • • • > ftw) • (A - 9) 

where the sum is taken over all possible permutations n of the set {1, . . . ,«}. These equations enable us to calculate the density 
contrast in the 71-th order of perturbation theory by using the iterative equations for the coupling functions. 
The calculation of the second-order coupling functions is straightforward. The result is 

Us), , 5 2 (q r q 2 ) 2 l q l -q 2 (qi q 2 \ (s) 3 4 ( gl ■ q 2 ) 2 l q l -q 2 (q 1 q 2 

F 2 (ft,ft) = = + = — 2-^— + ~ — + — , G 2 '(q u q 2 ) =- + + — + — 

7 7 2 qiq 2 \q 2 q\ 7 7 2 4142 \# 2 <7i 



The third-order coupling function is given by 

Fliquqi, qi) = 7§{ 7a (ft> 923)^2(? 2 . ft) + 2 ^(?1'923)G2(?2> ft) + I7a(?12> ft) + 2/3(? 12 , ft)]G 2 (0, , ft)} . ( A - 1 1) 

where = ^ ( + ^ ■. Employing Eq. ( | A.9| >, we find the symmetric function 



^(ft, ft, ft) = ^Kft, ft 3 )^ 2 S) (ft, ?3) + »(ft> 9n) F f(9v93) + <*(ft, q\ 2 )Ff{qi,q 2 )} 

+ |^(ft , ft 3 )G 2 s) (ft, ft) + Aft, ft 3 )G 2 s) (ft, ft) + Aft. ft 2 )G 2 s) (ft , ft)] 
7 

+ g^Kfta- ft) G 2 (?1> ft) + a ^13' ft) G 2 Aft'ft) + «(ft 3 , ft^fCft, ft)] ■ (A. 12) 

From now on the symmetry superscript "(s)" will be omitted because we will only deal with symmetric coupling functions. For the 
calculation of the trispectrum in the halo model approach as described in Sect. 4.5 one needs perturbation theory. More precisely, 
we need the subsequent components 

7 

^i(ft, -ft, ft) = ^Wft, ft) F 2(-ft, ft) + a(-q l ,q + )F 2 (q 1 ,q 2 )] 
4 

+ ^[^(ft, ft)G 2 (-ft, ft) +j8(-ft, ft)G 2 (ft, ft)] 
7 

+ 54^-' ft) G 2(-ft>ft) + -ft)G 2 (ft, ft)] , (A.13) 

where we have defined the difference vector q = q 2 - q { and the sum of the vectors q + = q 1 + q 2 . 

We already mentioned in the previous section that it is possible to find a solution for an arbitrary cosmology if one makes a small 
approximation. In the literature one can find closed solutions for the second- and third-order coupling functions. The second-order 
coupling function changes to 

F 2 (ft,ft) = i(l + e)4^f^ + ^ + f|-f)^#, (A.14) 
2 2 qiq 2 \q 2 q x j \2 2) qfq* 

where e * (3/7)Q m 2/63 for £l m > 0.1 |Bernardeau et al.|2002| l. For our fiducial choice of fi m = 0.3, we get Q m 2/53 ~ 1.039. Thus, 
within a few percent correction to the first and last term, the coupling functions are independent of cosmological parameters. 
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A.2. Correlation functions 
A.2.1. Bispectrum 

The dark matter bispectrum is defined as 

<5(Jfci)3(*2)S(* 3 )>c = (2jr) 3 6 D (k l2 3)B(k u k 2 ,k 3 ). (A.15) 

Since the connected bispectrum vanishes for Gaussian random fields, it is the first intrinsically non-linear moment. Inserting the 
perturbative expansion ( |A.4[ > for each term results generally in an infinitely large sequence of correlators. The lowest non-vanishing 
order is the so-called tree-level contribution to the bispectrum. We find for the correlator in tree level 

^(fcOfeWXree = (kfotfltetflte)) + 01 (*3)> + <?1 (*I )*1 (* 2 )§2(* 3 )> ■ (A. 16) 

Replacing the second-order density contrast with Eq. |a3J results in 

<5 2 (*i)5i(* 2 )5i(*3)> = J Jd 3 92 fe(*:i- 9l - 92 )F 2 ( ?1 , 92 X§i( 9l )§i(?2)^i(*2)§i(*3)> 

= (2tt) 3 5 d (A: 12 3) 2F 2 (A: 2 , k 3 )PUk 2 )PUk 3 ) , (A. 17) 

where we applied Wick's theorem to express the four-point correlator of Gaussian fields in terms of products of power spectra, and 
performed the two integrations over the Dirac delta distributions. The results for the other two terms of the tree-level bispectrum are 
simply obtained by permutations of the arguments. Finally, the tree-level bispectrum is given by 

B vt (k u k 2 ,k 3 ) = 2F 2 (k u k 2 )P l P 2 +2F 2 (k u k 3 )P 1 P 3 +2F 2 (k 2 ,k 3 )P 2 P 3 . (A.18) 

The factor 2 follows from using the symmetrized version of the second-order coupling function. 

A.2.2. Trispectrum 

The dark matter trispectrum is defined as the connected four-point function in Fourier space: 

<5(*i )~5(k 2 )6(k 3 )6{k 4 )) c = (27ifd D (k l234 )T(k u k 2 ,k 3 ,k 4 ), (A. 19) 

where £1234 = k\ + k 2 + k 3 + k 4 . We find that there are two different non-vanishing contributions to the tree level: 

(S(k\)d(k 2 )d(k 3 )S(k4)) tree = <§ 2 (*i)§ 2 (* 2 )§i(*3)§i(*4)> + . . . (6 terms) 

+ <§3(fci)§i(fc 2 )^fc3)§i(fct)> + • . . (4 terms) . (A.20) 

In total we find 6 terms of the first type and 4 terms for the second type, where the rest is obtained by permutations. All other 
contributions either vanish or are built of higher-order terms. Note that for the second type of terms we need the results from 
perturbation theory up to the third order. The calculation of each term is a tedious but straightforward calculation. We obtain for the 
first term of the expansion 

<5 2 (A:i)§2(*2)§i(*3)§i(*4)> = (2n) 3 6 D (k 1234 )4P 3 P 4 [P 13 F 2 (k 3 , -k 13 )F 2 (k 4 , -k 24 ) + P 14 F 2 (k 3 , -k 23 )F 2 (k 4 , -k 14 )] , (A.21) 

where the six-point correlator resolves into 15 terms consisting of power spectra products. Performing the integrations over the 
arising delta functions yields in the end 8 different terms. Similarly, we find for the second type of terms 

(8 3 (ki)Si(k 2 )Si(k 3 )di(k4,)) = (2n) 3 6 D (k l234 ) 6F 3 (k 2 ,k 3 ,k 4 )P 2 P 3 P 4 . (A.22) 

The other terms are easily obtained by permutations, however, we present here the complete result to avoid confusion with a 
shorthand notation that is introduced afterwards. The trispectrum of cold dark matter is in first non-vanishing order given by (Fry 
[1984) 1: 

r pt = 4T a + 6T b , (A.23) 

where 

T a = P1P2 [PnF 2 (k u -k l3 )F 2 {k 2 , k l3 ) + P 14 F 2 (k u -k 14 )F 2 (k 2 , k u )] 
+ P1P1 [PnF 2 (ku-k 12 )F 2 (k 3 , k l2 ) + Pi 4 F 2 (k u -k u )F 2 (k 3 , *„)] 
+ PiP 4 [Pi 2 F 2 (k u -k 12 )F 2 (k 4 , k l2 ) + P 13 F 2 (k u -k 13 )F 2 (k 4 , k 13 )] 
+ P2P3 [P2iF 2 (k 2 , -k 2l )F 2 (k 3 , Jfcai) + P 24 F 2 (k 2 , -k 24 )F 2 (k 3 , k 24 )] 
+ P2P4 [P2iF 2 (k 2 , -k 21 )F 2 (k 4 , fci) + P2 3 F 2 (k 2 , -k 23 )F 2 (k 4 , k 23 )] 

+ P 3 P 4 [P 3l F 2 (k 3 , -k 3l )F 2 (k 4 , k 3l ) + P 32 F 2 (k 3 , -k 32 )F 2 (k 4 , k 32 )] , (A.24) 



and 



T h = F 3 (k u k 2 , k 3 )P x P 2 P 3 + F 3 (k 2 , k 3 , k 4 )P 2 P 3 P 4 + F 3 (k 3 , k 4 , k{)P 3 P A P x + F 3 (k 4 , k x ,k 2 )P 4 P x P 2 , (A.25) 
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where P, = P^ik), Pij = PimQki + kj\) and k u = £, + kj. 

For the covariance matrix one only needs the parallelogram configuration. This imposes the condition k 2 - -k\ and k\ = -£3 
on the wave-vectors. In this case Eq. ( |A.23) > simplifies to 

7p t =4P 2 1 {[F 2 (k u -k + )fP + + [F 2 (k u k-)] 2 P-}+4Pl{[F 2 (k 3 ,-k + )?P + + [F 2 (k 3 ,-k_)] 2 P-} 
+ 8P,P 3 [F 2 (k u -k + )F 2 (k 3 , -k + )P + + F 2 (k u k-)F 2 (k 3 , -k-)P.] 

+ l2[p 2 1 P 3 F 3 (ku-k 1 ,k 3 ) + P 1 PlF 3 (kuk 3 ,-k 3 )\ , (A.26) 

where k = k 3 - k\, k + = k\ + k 3 , P = P^Qk-l) and P + = P\ m {\k + \). Consequently, we define the non-Gaussian contribution to 
the covariance in tree-level perturbation theory as 

C pt = C pt (4 £j) = - —— —— dw z-T pt \—, , — , ,w\ , (A.27) 

A J|;,|Ef, A r (£i) J\i 2 \ e ( i A x {(j) (2n) 2 J w 2 \w w w w ) 

where A denotes the survey area, A r (Z) the integration area and G(w) the lensing weight function (see Eq. |7J. More details on the 
notation can be found in Sect. [3] 
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